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ABSTRACT 

This paper examines how type I planet migration is affected by the presence of 
turbulent density fluctuations in the circumstellar disk. For type I migration, the 
planet does not clear a gap in the disk and its secular motion is driven by torques 
generated by the wakes it creates in the surrounding disk fluid. MHD turbulence 
creates additional density perturbations that gravitationally interact with the 
planet and can dominate the torques produced by the migration mechanism itself. 
This paper shows that conventional type I migration can be readily overwhelmed 
by turbulent perturbations and hence the usual description of type I migration 
should be modified in locations where the magnetorotational instability is active. 
In general, the migrating planet does not follow a smooth inward trend, but 
rather exhibits a random walk through phase space. Our main conclusion is 
that MHD turbulence will alter the time scales for type I planet migration and 
- because of chaos - requires the time scales to be described by a distribution of 
values. 

Subject headings: planets: formation - extrasolar planets 

1. Introduction 

Ever since their initial discovery (Mayor & Queloz 1995, Marcy & Butler 1996), 
extrasolar planets have been found in unexpectedly short-period orbits. Since theories of 
planet formation suggest that planets must form further out in the disk (e.g., Lissauer 
1993), these findings strongly indicate the necessity of planetary migration (e.g., Lin, 
Bodenheimer, & Richardson 1996). A variety of migration scenarios have been suggested. 
Type I migration (Ward 1997) occurs when a planet (or protoplanet) is embedded in a 



- 3- 



circumstellar gaseous disk and the planet's mass is too small to clear a gap in the disk. In 
this case, the planet drives wakes into the background fluid of the disk and these wakes, 
in turn, exert torques on the planet. These torques are generally greater in the inward 
direction than in the outward direction, and a net inward migration occurs. 

These same disks are expected to contain magnetic fields that are subject to a robust 
instability (Balbus & Hawley 1991), for regimes in which an ideal-MHD description is 
adequate. The nonlinear outcome of this magneto-rotational instability is MHD-driven 
turbulence. Indeed, MHD-driven turbulence is commonly invoked as a primary source of 
disk viscosity which can drive mass accretion through the disk. The usual scenario for type I 
migration, however, implicitly assumes a smooth background, uninterrupted by turbulence. 
Given that such turbulence, along with associated stochastic density fluctuations in the 
disk gas may often be present, we re-examine the type I migration scenario by considering 
the dynamical evolution of embedded protoplanets that are torqued by transient density 
perturbations arising from MHD turbulence. 

The effects of large-scale turbulence on type II migration, where the planet clears a gap 
in the disk (e.g., Goldreich & Tremaine 1980; Lin & Papaloizou 1993), has been explored 
recently (Nelson & Papaloizou 2003; Papaloizou & Nelson 2003; Winters, Balbus, & Hawley 
2003). In a related vein, Terquem (2003) has considered the possibility of stopping inward 
migration by a strong toroidal magnetic field. These previous papers show how turbulence 
and magnetic fields can produce important modifications to the standard migration picture. 
This present work carries out a complementary analysis for type I migration. 

Our approach to the problem is as follows. We first compute a representative 
three-dimensional MHD simulation (§2) in a circumstellar disk and use the results to 
quantify the spectrum of nonlinear density perturbations that are produced by MHD-driven 
turbulence. These density perturbations provide a stochastic time-varying gravitational 
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torque on any embedded protoplanets that are orbiting in the disk. Our next step is to 
construct a parametric model which characterizes the essential aspects of the spectrum of 
turbulent density fluctuations that are observed in the three-dimensional simulation. This 
model is described in terms of a time- varying non-axisymmetric potential function(§3). 
The resulting prescription for the perturbations is then applied to two-dimensional 
hydrodynamical simulations of circumstellar disks to study type I protoplanetary migration 
(§4). With this scheme, we can examine the effects of the turbulent perturbations on 
the migration of bodies that are not sufficiently massive to open gaps in a protostellar 
disk, and we can obtain an overall initial survey of the Type I migration problem. Our 
approach allows for a rapid survey of the important regions of parameter space, and is 
complementary to full-scale 3D MHD simulations which trace the evolution of a planet in 
an ionized circumstellar disk. Because the saturation of the magnetorotational instability 
occurs at nonlinear amplitudes, MHD-driven turbulence naturally produces large density 
fluctuations, which in turn produce torques that generally dominate those produced by the 
planetary wakes of objects with M < 10M e . We quantify this threshold by deriving an 
analytic estimate of the fluctuation amplitudes that are required for type I migration to be 
highly distorted (§5) and then discuss the longer term behavior (§6). 



2. MHD Simulations 

Our first step is to conduct a representative three-dimensional simulation to evaluate 
the expected properties of MHD turbulence in a model cicumstellar disk. We use the 
"Nirvana" MHD code (Ziegler & Riidiger 2000) as modified by Steinacker & Papaloizou 
(2002). We work in cylindrical coordinates (z,R,(j>), and use a system of units in which 
the gravitational constant, G, and the central stellar mass, M, are both unity. Our disk 
model has a nearly Keplerian radial domain that extends from an inner boundary at R—1.5 
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out to R=3.5, and covers an azimuthal slice that subtends n/3 radians. The equation 
of state is taken to be locally isothermal, with P(R) = p(R)c s (R) 2 (see e.g. Laughlin 
& Bodenheimer 1994). This approximation is reasonable for optically thick disks whose 
local radiative diffusion time is shorter than the local dynamical (orbital) timescale. The 
gravitational potential from the central star is specified to depend only on the cylindrical 
radial coordinate $ = —GM/R. The model therefore does not consider the vertical density 
gradient that results from the z component of the gravitational force. This approximation 
simplifies the physical situation, allowing the use of a modest number of zones in the 
vertical direction. An obvious future refinement will be to study vertically stratified models. 
The reference simulation described here uses 334 radial zones, 72 azimuthal zones, and 40 
vertical zones. 

In the nearly Keplerian region of our initial equilibrium model, the disk density varies 
inversely with radius, p(R) = p (R /R), the soundspeed is given by c s (R) = C ^/Rq/ R, 
and the rotational velocity is v<p(R) = ^/GM/R — 2Cq^R /R. For numerical reasons, the 
disk also contains inner and outer radial boundary domains as described in Steinacker & 
Papaloizou (2002). 

We adopt Cq = 0.1, along with a subthermal poloidal seed field having zero 
net flux through the disk. With these starting parameters, the disk experiences the 
magnetorotational instability, progresses through the channel phase, and develops 
turbulence as described in detail by Steinacker & Papaloizou (2002). This configuration of 
initial conditions was found by Steinacker and Papaloizou (2002) and Hawley (2001) (among 
others) to lead to the smallest overall amplitude in the turbulent density fluctuations. Initial 
magnetic configurations with net flux threading the disk can lead to turbulent fluctuations 
that are up to two orders of magnitude larger. This turbulence persists for the duration of 
the simulation, which was followed to 200 orbits at the outer disk edge. 
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Figure 1 shows a map of the surface density distribution cr(R, <fi) in the disk after the 
turbulence has entered its quasi-steady fully developed phase. The order-unity surface 
density fluctuations that are present in the disk provide the stochastic gravitational torques 
the are capable of severely modifying the type I migration process. The corresponding 
global Fourier amplitudes 

Jo J Rin e{r,<t>)rdrd<t) 

of the disk fluctuations shown in Figure 1 are plotted as the filled black circles in Figure 
2. The perturbation amplitude is strongest for the lowest azimuthal wavenumber (m = 6) 
that fits in the n/3 azimuthal domain. As the azimuthal mode number m = 6n is increased 
through n=2,3,4, and 5, the Fourier amplitudes show a gradual decline. The concentration 
of power in lower-order Fourier modes is a characteristic of MHD-generated turbulence, and 
has been previously observed by other authors, including Armitage (1998), and Hawley et 
al (1995). 



3. Parameterization of the Turbulence 

In order to survey the effects of MHD turbulence on type I planet migration, we 
have chosen to construct a simple heuristic model to describe the perturbations. By 
characterizing the turbulence arising in the full MHD treatment in a simple manner, we can 
use the resulting prescription to specify the development of time-varying disk fluctuations 
in two-dimensional planet migration simulations. The alternative (and more rigorous) 
approach to the problem would involve running an expensive three-dimensional MHD 
simulation with every planet migration calculation. Each long-duration run of the MHD 
code requires, for example, 30,000 CPU hours on an Origin 3800 parallel processing facility 
(Nelson & Papaloizou 2003). Running enough simulations to sample the parameter space 
of interest is thus overly time-consuming. 
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We are primarily interested in the gravitational forces exerted on a migrating planet by 
the density perturbations arising from the MHD instabilities. We thus need a description of 
the potential perturbations <E>, which can be differenced to compute forces that act on the 
fluid elements; the gravitational forces from the resulting density enhancements then act on 
the planet itself. Specifically, we adopt a simple heuristic model for the potential <3> induced 
by the turbulent perturbations: 

$ = ~2 cos[m6 -ip- Q c t)\ sin[vr— ] . (2) 

With enough modes of this form, one can characterize the spectrum of density fluctuations 
produced by the MHD turbulence. Given the finite resolution of our codes, the computed 
MHD perturbations can be reasonably modeled by applying 50 modes at any given time. 
The modes come in and out of existence with the time dependence specified above, where 
t = t — t . An individual mode begins at time t an d has faded away after a time At 
= tf — to- After the mode is gone, a new mode is generated. With the proper choice of 
the distributions of the remaining parameters (see below), we can reasonably mimic the 
power spectrum of density fluctuations in the actual MHD turbulence computed by the full 
three-dimensional code (see Figure 2). 

The following parameter choices are found to work: The modes are centered at radial 
location r c and initial angular location (phase angle) if. The location of the mode r c is 
chosen from a random (uniform) distribution, as is the phase angle ip. The extent of the 
mode is determined by the azimuthal wave number m which is chosen to be distributed 
according to a log-random distribution for wavembers in the range 2 < m < N grid/8 
(where N gri d is the number of azimuthal zones in the simulation, typically 512 here). We 
find that the choice of a log-random distribution of m gives the best approximation to the 
density fluctuation spectrum produced by the MHD simulation; this spectrum of m provides 
power on all spatial scales. With m specified, the mode extends for a distance 2nr c / m along 
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the azimuthal direction. The radial extent is then specified by choosing a = itr c /Am so 
that the mode shapes have roughly a 4: 1 aspect ratio (this statement is not precise because 
the profiles differ in the radial and azimuthal directions). Notice also that the initial mode 
shapes are not tilted as would be expected from the usual spiral patterns seen in disks; 
here, the durations of the modes are generally much longer than the time required for the 
Keplerian shear to adjust the orientations, so we can rely on the shear to take over. 

The modes first appear at a time to, which is chosen so that the simulation contains 50 
active modes at all times. The pattern speed Vt c = v kep i er /r c in the time dependent factor 
allows the mode center to travel along with the Keplerian flow. The temporal duration of 
the mode At = tf — 1 is taken to be the sound crossing time of the mode along the angular 
direction, i.e., At = 2rrr c /(mas). 

Finally, the modes have an overall amplitude A and an individual amplitude £. As 
written, the amplitude A has units of [A] = £ 5 / 2 t~ 2 (where £ is a length scale and t is a time 
scale) so that the overall expression has units of potential. The dimensionless variable £ has 
a gaussian distribution with unit width, whereas the overall amplitude A has a fixed value 
for all of the modes. We can find the amplitude A that reproduces the level of perturbations 
found in the MHD simulation. Alternately, we can treat the amplitude A as an adjustable 
parameter and study the effects of these turbulent fluctuations on planet migration as a 
function of the amplitude (see below). 

Figure 2 shows the power spectra from both the full MHD treatment of the previous 
section and the parametric model developed in this section. The two spectra are in good 
agreement. As a result, we employ this approach to specify the perturbations for migration 
calculations, as presented in the following section. 
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4. Migration with Turbulent Perturbations 

With the perturbations due to MHD turbulence specified (as described in §3), we 
have performed fluid dynamic simulations of the migration of a terrestrial planet (or core 
of a giant planet) in the circumstellar disk. These simulations use a two-dimensional fluid 
dynamics code where we add the turbulent perturbations according to the prescription 
described above. 

Our two-dimensional code (Laughlin 1994) is based on second-order van Leer advection 
(e.g. Stone & Nornam 1992). We generally use a grid of 256 evenly spaced radial zones 
and 512 azimuthal zones. The initial model has the same surface density and sounspeed 
distributions and the same active radial domain as was employed in the three-dimensional 
simulation. We minimize reflection by implementing sponge boundary conditions at the 
radial edges of the computational domain. In the = 12 radial zones interior to £ = 1.6, 
and in the rid = 12 radial zones exterior to £ = 3.4, we smoothly impose an admixture 
of the initial equilibrium solution into the hydrodynamically computed variables x hydro (£). 
That is, at the inner edge we have, 

X (Zj) = (- ^Kydrofe) + (— jKquilfe) • (3) 

This damping is applied at a cadence tdi = A£(ndi)/cg at the inner edge, and 

t do = A£(n<2 ) / c g at the outer edge, where A£(j) is the radial width of the uniformly spaced 

zones. 

The planet mass is taken to be either 1 or 10 M e , and is initially placed in a circular 
orbit of semi-major axis (radius) a = 2.5 AU. (Note that all of these radii can be scaled). 
Over this restricted annulus covered by the hydrodynamical domain, the initial disk mass 
set to be a small fraction of the stellar mass, i.e., = 0.02. The subsequent evolution 

of a 1M® planet in the a — e plane is shown in Figure 3 for three different amplitudes of 
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the turbulent perturbations: A = 5.0 x 10~ 6 , A = 5.0 x 10~ 5 , and A = 5.0 x 10~ 4 . The 
corresponding surface density perturbations are depicted in Figures 4, 5, and 6. 

For small fluctuation amplitudes A, type 1 migration proceeds in largely unperturbed 
fashion. The planet slowly and steadily spirals inward. The fluid dynamic simulations only 
follow the planet for ~ 100 orbits in this class of simulations (see Figs. 3a and 4a). When, 
on the other hand, the fluctuation amplitude is large, the chaotic perturbations produced 
indirectly through turbulence dominate the torques produced by steady type I migration. 
In this regime, the evolution of the planet is highly chaotic (see Figs. 3c and 4c). In between 
these extremes, there exists a critical fluctuation amplitude, the value of A at which type I 
migration changes its character, but a general inward migration trend is still evident (Figs. 
3b and 4b). The critical amplitude determined empirically from these two-dimensional disk 
simulations is Ae ~ 1 x 10 -5 . In the following section, we derive an analytic estimate for 
this critical value. 



5. The Critical Amplitude 

Equation (2) gives the potential produced by the turbulent fluctuations. The gradient 
of this potential provides a forcing term that acts on the gas in the disk (but does not 
act directly on the planet) and thereby produces corresponding fluctuations in the surface 
density of the disk. To obtain a rough understanding of the critical fluctuation amplitude 
A c required for turbulence to overwhelm the standard type I migration torques, we assume 
here that the gravitational potential of the perturbed gaseous disk (where this potential 
does act on the planet) has the same form as that of equation (2). However, the amplitude 
of the potential acting on the planet will, in general, be reduced from that acting on the 
gas by a factor T < 1 (which we estimate below). The torques exerted on the planet by the 
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surface density perturbations in the disk will thus have the form 



mAtTMp 



^ = s--^ e -(r-r c )V^ gin [ m ^ _ ^ _ sin [ n i/At] , (4) 



r 



1/2 



where Mp is the mass of the planet and all the other quantities are defined previously. At 
any given time, the disk will contain 50 modes that exert torques of this form. To evaluate 
the effectiveness of these torques, we need to define an average torque strength. First, we 
integrate over time to obtain the net change in angular momentum per mode, i.e., 

A J = r dt r = mA ^ /2 Mp e-^f^ * sin^-y^-siny. 
J r 1 ' 2 7r a 2 — 1 

where the parameter a is defined to be a = (mVL — VL c )At/rc. The mean torque per mode t±, 

averaged over the lifetime of the mode, is thus given by 

A J ATMp mi _ {r _ rc) 2,2 sin(a7r - tp) - siny? 

n = — = — -^T where T = — e {r Tc> /<T , (6) 

At r v i l 7r a 1 — 1 

where the variables r c , ip, ^, and m are chosen randomly as described in §3. 

At any given time, the disk contains 50 modes and hence 50 values of the function T . 
Since these torques can either be positive or negative, the 50 modes will tend to cancel each 
other out, and the net effect forcing effect will be relatively small. We thus need to compute 
the forcing function JF 50 defined by 

50 

•^50 = ^ Fj , (7) 

J'=l 

where the Tj are sampled in the same way that the numerical simulations sample parameter 

space. The function JF 50 should average to zero (or some low value) because the torques 

can be positive or negative. However, the RMS value of JF 50 , denoted here as (.F50), should 

provide a good measure for the effective strength of this torque. Numerically sampling the 

function shows that (JF 50 ) ps 0.29. Putting these results together, we find that the average 

torque Tp exerted on the planet through the action of turbulence is given by 

, , ATMp ATMp 
rT=(^o)- TJ f^0.29- TJ f. (8) 
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To estimate the response amplitude T of the disk (due to the potential induced by 

turbulent fluctuations), we use a simple linearized WKB treatment of a two-dimensional 

disk (Shu 1992). Unlike the usual spiral density wave theory, however, we include only the 

forcing due to the fluctuation potential $ and neglect the self-gravity of the disk. In the 

WKB limit, we can solve for the radial part of the reduced surface density perturbation S 

in terms of the forcing potential $ to obtain 

= k 2 * d $ 

{u - mtt) 2 -k 2 - k 2 a 2 s ' 1 ' 

where the epicylic frequency k = Q for a Keplerian disk. Next, we need to use the 

relationship between the surface density perturbation S and the corresponding gravitational 

potential V. In the WKB limit, this relation has been calculated previously (e.g., Shu 1992) 

and can be written in the form 

Equating the expressions for the surface density response to the turbulent fluctuations (eq. 
[9]) and the corresponding perturbation in the gravitational potential (eq. [10]), we can 
solve for the response factor T, i.e., 

r ~ V ~ 2irGo~dk 

~ $ ~ (u - mny -n 2 - k 2 a 2 ' y } 

To evaluate this expression, we must make further simplifying assumptions. The turbulent 
fluctuations induce modes with relatively large wavenumbers m = 2 - 32, so we can 
approximate the denominator as \V\ ~ m 2 Vl 2 = m 2 GM*/r 3 . Given the radial dependence of 
the turbulent fluctuations (eq. [2]), we can estimate the radial wavenumber k in the WKB 
limit: kr = (r/&)(d&/dr) ~ 2r 2 /a 2 ~ 32m 2 /n 2 . The expression for T thus becomes 

Tio A r 2 64 



M* 7T 2 



0.06 . (12) 



Next, we need to compare the average torque produced through turbulent fluctuations 
with that produced by the usual type I migration mechanism (Ward 1997). This torque 77 
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can be written 

n = Pr(M P /M*) 2 (nay)(rn) 2 (r/H) 3 , (13) 

where H is the disk scale height, a A is the disk surface density, and f3j is a dimensionless 
amplitude. Previous calculations suggest that (3j ~ 1CT 2 (Ward 1997). 

By equating t t and tj, we can solve for the critical amplitude A c for which the torques 
due to turbulence become larger than those of standard type I migration. We find 

A c = -^(3j(M P /M^(r/H) 3 r^ 2 (rn) 2 » 1 x 10" 5 , (14) 

where the numerical value is obtained using the same parameters as in the numerical 
simulations presented here (for a one Earth mass planet). This analytic estimate for the 
critical amplitude is in good agreement with that observed in the numerical simulations 
themselves - which show an empirically determined value of A E ~ 1 x 1CT 5 (see §4). 

Before leaving this section, we note that the torque produced by type I migration (eq. 
[13]) and that produced indirectly through turbulence (eq. [8]) scale differently with the 
mass of the planet. Small planetesimals, with masses m <C M E , will be tossed around 
violently because of MHD turbulence. On the other hand, sufficiently large planets can 
be impervious to its effects (for sufficiently small amplitudes A). If we assume that the 
maximum fluctuation amplitude A max r* r 1 / 2 (rfi) 2 , then we can define a critical mass scale 
Mcp, such that larger masses will always be dominated by type I migration torques. This 
mass scale is given by M C p/M* = (64/V 2 /3j) (r/H) 3 r* 0.04. Secondary bodies with such 
large masses are generally brown dwarfs and not planets. Furthermore, such large bodies 
tend to open up gaps and experience type II migration (Goldreich & Tremaine 1980; Lin & 
Papaloizou 1993; Ward 1997). As a result, there is no regime of parameter space for which 
type I migration torques necessarily dominate those produced by MHD turbulence. 
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6. Long Term Evolution 

Given our analytic description of the effects of turbulence, we can understand the 
possible range of migration behavior over longer time intervals. This analytic approach 
is necessary, as it is difficult for numerical simulations to remain stable over millions of 
orbits, the dynamical time frames of the migration process. As we discuss below, several 
complications arise in the long term. 

First, even if the amplitude of the turbulent fluctuations is that required for the 
resulting torques to have the same amplitude as those of type I migration (i.e., if A = Ac), 
the two sources of torque have different accumulative effects. The type I migration torque 
Tj is a secular torque that consistently moves the planet in one direction. The net turbulent 
torque tt can be either positive or negative and changes its value over a sampling time that 
we denote here at t smp . The changes in angular momentum accumulate in random walk 
fashion. The change in angular momentum over the integrated time interval t = Nt smp , for 
the two contributions, can thus be written in the form 

(AJ)j = Nnt smp and (AJ) T = ^r T t smp . (15) 

The sampling time t smp is the time required for the disk to produce an independent 
realization of the turbulent torques. The lifetime of a mode with wavenumber m is At 
= 2nr I '(mas), which can be written in the form At = (2n/Q)(r/mH). The first factor 
is the orbital period P orb . The second factor is of order unity: The disk generally has 
r/H pa 10 and the wavenumber m lies in the range 2 < m < 32. Since the wavenumbers 
m are logarithmically spaced and the lowest m modes live the longest (and are near the 
peak of the power spectrum - see Fig. 2), the appropriate value of m to use in estimating 
the sampling time lies near the low end of the range. As a result, we expect £ smp ~ P OT b- 
In order for the turbulent fluctuations to dominate over a specified time interval tf, the 
amplitude of the turbulent fluctuations must be larger than the critical value A c by an 
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additional factor, i.e., 

A C . (16) 

A sufficiently large amplitude can keep the turbulent fluctuations dominant over the 
entire lifetime of the disk. The expected disk lifetime tdisk (specifically, the time over 
which planet migration takes place) is of order 1-3 Myr (Haisch, Lada, & Lada 2001), 
while expected timescales for Type I migration of individual (Earth mass) planetesimals 
are considerably shorter, of order 10 5 years (Ward 1997). The orbit time P D rb is of order 
1 -3 yr. So this problem contains another critical amplitude, namely that required for 
the turbulent fluctuations to overwhelm type I migration torques over a 10 5 year secular 
migration timescale. This second critical amplitude A 2 is estimated to be 

A2 ~ C-^-Y^Ac - 0.003. (17) 

^-P orb 7 

A useful benchmark is the number of steps N (or, equivalently, the time scale iVP or b) 
required for the random walk process to completely change the angular momentum. In 
other words, we set (AJ) T = \fGM~r- in equation (16) and solve for N to find 



N 



,5/2^2 ; 



-2 



(18) 



Even for fluctuations with the critical amplitude A = A c , this number is rather large: 
N ~ 3 x 10 8 , with a corresponding time scale of 300 Myr. On the other hand, for A = A 2 , 
the required number of steps falls to iV ~ 3 x 10 3 , with a time scale of ~ 3000 yr. 

Another complication is the radial dependence of the critical amplitude Ac- As shown 
by equation (14), this ratio scales as A c oc r 3 ~ 3<? / 2 , where q is the power-law index of the 
disk temperature profile. Although the temperature will not be purely a power-law, q is 
nonetheless expected to lie the range 1/2 < q < 3/4, so that the index 3 — 3g/2 = 15/8 — 
9/4^2. As a result, the amplitude of the turbulent modes required to dominate the type 
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I migration torques is a rapidly increasing function of radius r (decreasing as the planet 
moves in). As the planet migrates inward, A c oc r 2 decreases and turbulence tends to 
dominate to an ever greater extent. On the other hand, if MHD turbulence shuts down in 
the inner disk due to magnetic decoupling (e.g., Fromang, Terquem, & Balbus 2002), then 
standard type I migration is resumed. 

7. Discussion and Conclusion 

In this paper, we have generalized the scenario of type I planet migration to include the 
effects of fluctuations produced by MHD turbulence. To characterize the turbulence itself, 
we have run three-dimensional MHD simulations to study the onset and character of the 
turbulence, and used the results to specify the power spectrum of turbulent fluctuations. 
These perturbations were then parameterized and used in a second set of two-dimensional 
numerical simulations that study the migration phase itself. Finally, we have developed an 
analytic description of the turbulent fluctuations, the torques they produce, and their long 
term effects on the migration process. Our specific results are summarized below: 

MHD turbulence produces a full spectrum of fluctuations that peaks at relatively low 
azimuthal wavenumbers m ~ 1 — 6 (Figs. 1 and 2). These fluctuations can be characerized 
by potential fluctuations of the form given by equation (2). The formulation developed here 
may be useful in many other astronomical applications because it allows the fluctuations 
produced by turbulence to be effectively modeled with a minimum of computational effort, 
and even allows for many results to be determined analytically (see §5). 

The effects of turbulence on type I migration depend sensitively on the amplitude of 
the turbulent potential perturbations. Further, we have found the critical threshold Ac 
for this amplitude, using both two-dimensional numerical simulations of the migration 
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process (Figs. 3 and 4) and an analytic treatment (eq. [14]). For perturbation amplitudes 
below this threshold (A < A c ) type I migration proceeds smoothly inward, as envisioned 
by the original scenario. For small but finite amplitudes, the turbulence leads to random 
fluctuations superimposed on the smooth inward migration, but the general trend remains. 
If the amplitude exceeds the threshold (A > A c ) , the perturbations due to turbulence 
produce torques that are stronger than those induced by the planetary wake, and the 
evolution changes dramatically. Instead of displaying a smooth, inward progression, the 
planet exhibits random walk behavior. The planet can even move outward when the 
turbulence amplitude is large enough. 

The critical amplitudes are A c ~ 1 x 1CT 5 for a 1 Earth mass planet and A c ~ 1 x 10~ 4 
for a 10 Earth mass planet. These critical amplitudes are found both in the two-dimensional 
numerical simulations of type I migration and by the analytic treatment of §5. For 
comparison, the expected amplitudes for the perturbations are A ~ Z ~ 1CT 3 , comfortably 
larger than the critical amplitudes needed to change the migration scenario. As a result, 
we expect MHD turbulence to change the character of the type I migration process under 
typical disk conditions. 

The main conclusion of this paper is that the typical migration torques in a circumstellar 
disk will often be dominated by the perturbations due to MHD turbulence rather than the 
steady planetary wakes that drive type I migration. As a result, type I migration will often 
display far richer behavior than has been considered previously. When MHD instabilities 
are robust, migration will generally proceed in a highly chaotic fashion, with an element of 
a random walk behavior, although the underlying conventional type I secular torque acts to 
drive the planet inwards over the long term. The complications introduced by turbulence 
have two important implications: (1) The time scale for planet migration can be quite 
different, either longer or shorter, when turbulence is present. (2) The result of any given 
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migration episode in a circumstellar disk cannot be described by a single outcome; instead, 
due to chaos and extreme sensitivity to initial conditions, the result of a migration episode 
must be described in terms of a distribution of possible outcomes. In particular, the time 
scale required for planets to migrate inward will display a full distribution of values. 

This chaotic element - including the necessity of finding the full distribution of 
possible outcomes - has recently been emphasized for migration scenarios that involve 
planet-scattering (Adams & Laughlin 2003). Other recent work (Nelson & Papaloizou 
2003; Papaloizou & Nelson 2003; Winters, Balbus, & Hawley 2003) outlines the effects of 
turbulence and chaos for type II migration. Since this paper shows that type I migration 
can often be dominated by turbulence/chaos as well, nearly all possible mechanisms for 
planet migration are expected to display chaotic behavior. The ubiquity of chaotic processes 
during the formation and early evolution of planetary systems certainly contributes to the 
startling diversity that is observed among the known extrasolar planets 
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Fig. 1. — Surface density distribution in a 7r/3 wedge of a circumstellar disk resulting from a 
three-dimensional simulation of MHD turbulence. Gravitational torques arising from these 
fluctuations will stochastically alter the orbits of protoplanets embedded in the disk. 

Fig. 2. — Non-axisymmetric modal amplitudes observed in the numerical simulations. The 
solid black circles show the global fourier amplitudes C 6m measured in the MHD-turbulent 
three-dimensional simulation at the time plotted in Figure 1. The solid red triangles show 
the amplitudes C m observed in a two-dimensional hydrodynamical simulation in which a 
lM e protoplanet is embedded in a 0.02 M & disk, and perturbed as described by equation 2 
with A = 5.0 x 10~ 4 . The open red squares show the amplitudes C m observed in a similar 
simulation with A = 5.0 x 10~ 5 . The solid blue circles show the amplitudes C m observed in 
a simulation containing a 1M® protoplanet embedded in a disk with A = 0. 

Fig. 3. — Type I migration simulations using a two-dimensional fluid code and the para- 
metric treatment of the perturbations due to MHD turbulence. The three panels show time 
evolution of the migrating planet in the a — e plane for fluctuations with low amplitude (a) 
where A = 5.0 x 10~ 6 ~ Ac/2, threshold amplitude (b) where A = 5.0 x 10~ 5 ~ 5A C , and 
high amplitude (c) where A = 5.0 x 10~ 4 ~ 50A C - All of these simulations use a planet with 
one Earth mass. Planets with higher mass behave similarly if the amplitude of the turbulent 
fluctuations is increased accordingly (see text). 

Fig. 4. — Type I migration simulation using a two-dimensional fluid code and the parametric 
treatment of the perturbations due to MHD turbulence. Perturbations in the surface density 
of the disk are shown after 50 orbits of an embedded lM e protoplanet for the case of 
turbulent fluctuations with low amplitude A = 5.0 x 10~ 6 ~ A c /10. The color scale ranges 
from a(r, 0)/ < a(r) >=0.998 (black) to a(r, 0)/ < a(r) >=1. 00836 (orange). 
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Fig. 5. — Same as the previous figure, but with amplitude A = 5.0 x 10 5 ~ Ac- The color 
scale ranges from a(r,<j))/ < a(r) >=0.992 (black) to a(r,4>)/ < a(r) >=1. 00861 (orange). 



Fig. 6. — Same as the previous figure, but with amplitude A = 5.0 x 10 4 ~ 10Ac- The color 
scale ranges from a(r, 0)/ < a(r) >=0.899 (black) to a{r,4>)/ < a(r) >=1.0801 (orange). 
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